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Abstract. - Monte Carlo simulations of the two-dimensional XY model are performed in a 
square geometry with fixed boundary conditions. Using a conformal mapping it is very easy to 
deduce the exponent rj^ (T) of the order parameter correlation function at any temperature in the 
critical phase of the model. The temperature behaviour of ricr{T) is obtained numerically with a 
good accuracy up to the Kosterlitz-Thouless transition temperature. At very low temperatures, 
a good agreement is found with Berezinskii's harmonic approximation. Surprisingly, we show 
some evidence that there are no logarithmic corrections to the behaviour of the order parameter 
density profile (with symmetry breaking surface fields) at the Kosterlitz-Thouless transition 
temperature. 



The two-dimensional classical XY model is very famous in statistical physics, both for fun- 
damental reasons, as describing for example classical Coulomb gas or fluctuating surfaces, 
and because of its relevance for the description of real systems, e.g. Helium super-fluid films. 
The model undergoes a standard temperature-driven paramagnetic to ferromagnetic phase 
transition in d > 2, characterised e.g. by a power-law divergence of the correlation length 
near criticality, ^ ^ ' (^^e e.g. Ref. |l]) while it exhibits a rather different behaviour 
in two dimensions, where the existence of a low-temperature phase with conventional long 
range order is precluded by the Mermin- Wagner theorem Indeed there is no spontaneous 
magnetisation at any nonzero temperature in one or two-dimensional isotropic systems with 
short range interactions and having a continuous symmetry group, for example 0{n) spin 
models for n > 2, since order would otherwise be destroyed by spin wave excitations 
Such systems can however display a low-temperature phase with topological order and the 
defect-mediated transition is governed by unbinding of topological defects. At low temperature, 
the system is partially ordered, apart from the topological defects (vortices in the case of the 2d 
XY-mode\) which appear in pairs, in increasing number with increasing temperature. At the 
transition temperature, the pairs are broken and the system becomes completely disordered. 
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The low-temperature regime was investigated by Berezinskii Q in the harmonic approximation 
and the mechanism of unbinding of vortices was studied by Kosterhtz and Thoulcss ^ using 
approximate renormahzation group methods. For reviews, see e.g. Refs. 0, ^, 

This very pecuUar topological transition is characterised by essential singularities when 
approaching the critical point from the high temperature phase, ^ ~ e''«* , x ^ ^^^^ > while in 
the low-temperature phase, there is a line of critical points and the order parameter correlation 
function decays algebraically with an exponent fjaiT) which depends on the temperature. 
Consider a square lattice with two-components spin variables — {cos 9 j,^, sin 9 u,) with 
continuous 0(2) symmetry, located at the sites w of a lattice A of linear extent L, and 
interacting through the usual nearest-neighbour interaction 

- = ^'XlXl'^"' ■ = /i:^^cos(0^ ~9^+fi), (1) 

where /i labels the directions and /i is a unit vector in the /i— direction. The low-temperature 
or spin- wave behaviour is obtained in the harmonic approximation, after expanding the cosine, 
which should be justified at sufficiently low temperature where existence of order, at least at 
short range, is assumed: 



H Hn 



- J^EE(^"'- W- (2) 



w ^ 

Within this approximation, the two-point correlator becomes 

(0^, .4.) ^ Z-'Y[ I ^cos(0^, X e-*^S™S.''™-'™+'^^' 

W 

~ \WI-W2t^'^^'', (3) 

hence r](T) = T/2tt. The low-temperature phase thus has the characteristic features of a 
critical phase with local scale invariance, i.e. conformal invariance (invariance under rotation, 
translation and scale transformations, short-range interactions, isotropic scaling). Conformal 
invariance is known as a powerful framework for the description of two-dimensional critical 
systems . In the case of the XY model, the central charge is c = 1 in agreement with the 
continuous variation of critical exponents in the low-temperature phase. The scaling dimen- 
sions Xn,7yi of the primary operators On.m depend on a parameter p (called compactification 
radius). 



The spin operator is identified to the conformal operator Oi^ and the vortex operator to 
C'o.i ||ic| . This is coherent with the picture according to which in the high-temperature 
phase, the vortices are unbounded and produce the disordering of the system, the vortex 
chemical potential playing the role of a relevant operator. Then vortices and anti-vortices 
begin to bind, decreasing their relevance when the temperature decreases to become marginal 
at Tkt, where the condensation of defects occurs. There, xo,i must be 2 (for the RG 
eigenvalue y ~ d~ xqi to vanish), implying that p = 2. Hence at the KT transition, the 
exponent of the spin-spin correlation function decay takes the value predicted by Kosterlitz 
and Thouless, 77o-(7kT) = '^^1,0 = 1/4. The corresponding values in the critical phase below 
Tkt would in principle easily be deduced from equation (^), but the dependence of p on the 
temperature is not known. Unfortunately, not much is known in the intermediate regime 
between the spin wave approximation at low temperature and the Kosterlitz- Thouless results 
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at the topological transition and in fact the precise determination of the critical behaviour of 
the two-dimensional XY model in the critical phase remains a challenging problem. Many 
results were obtained using Monte Carlo simulations (see e.g. Refs. 0, |l|, ^ |l|, 0) 
at Tkt or high-temperature series expansions ]l7t , but the analysis was made difficult due 
to the existence of logarithmic corrections, e.g. x ^ ^^^^'^ (In^)^^'' in the high-temperature 
regime, or {(jiw^ ■ (jiw^) — \wi — u'2|~''''(ln \wi — W2\)^''^'^ at Tkt exactly and the value of r]a 
at Tkt was a bit controversial as shown in table 1 of referen ce p6[ . The resort to large-scale 
simulations was then needed in order to confirm this picture ]15|. 

In the following, we use a rather different approach which does not require so extensive 
simulations. Assuming that the low-temperature phase exhibits all the characteristics of 
a critical phase with conformal symmetry, we simply use the covariance law of n— point 
correlation functions under the mapping of a two-dimensional system confined inside a square 
onto the half-infinite plane. The scaling dimensions are then obtained through a simple 
power-law fit where the shape effects are encoded in the conformal mapping. This is the crucial 
point, since then even very small systems are well adapted to such fits, apart from irrelevant 
lattice effects. We note that mappings inside a square and various other geometries have been 
considered already, e.g. the moments of the magnetization and structure factors in 
the Ising model have been calculated in square systems. 

The order parameter correlation function in a square system, {4>wi ■ 4>W2)-, should in principle 
lead to the determination of critical exponents, but practically, it is not of great help since 
strong surface effects occur which modify the large-distance power-law behaviour. The corre- 
lation function is supposed to obey a scaling form which reproduces the expected power-law 
behaviour in the thermodynamic limit, for example {4>wi ■4>W2) = |''i'i~u'2|~''''/sq.(wi/L, W2/L), 
where /sq. encodes shape effects in a very complicated manner. Conformal invariance provides 
an efficient technique to avoid these shape effects, or at least, enable to include explicitly 
the shape dependence in the functional expression of the correlators through the conformal 
covariance transformation under a mapping w{z): 

(4i • <A-.> = \w'{z,)\--'\w'{z2)\-^'{$., ■ 4) (5) 

with 77^ — 2a; o- . In the semi-infinite geometry z = x + iy (the free surface being defined by the 
X axis) , the two-point correlator is fixed up to an unknown scaling function (apart from some 
asymptotic limits implied by scaling). Fixing one point zi close to the free surface (zi = i) of 
the half-infinite plane, and leaving the second point Z2 explore the rest of the geometry, the 
following behaviour is expected: (02^ • 4>z2)^oo ^ iui ^ 2/2)^^" V'(<^)j where the dependence on 
to — yiy2/\ zi — Z2 of the universal scaling function ip is constrained by the special conformal 
transformation [^0[ . In order to get a functional expression of the correlation function inside 
the square geometry w, one simply has to use the mapping w{z) which realizes the conformal 
transformation of the half-plane z = x + iy {0 < y < 00) inside a square w = u + iv of size 
Lx L < u < L/2, < v < L) with free boundary conditions along the four edges. This 

is realized by a Schwarz-Christoffel transformation |2l| 

w{z) = ^Fiz,k), z = sn(^^y (6) 

Here, F{z,k) is the elliptic integral of the first kind, sn (2Kw/L) the Jacobian elliptic sine, 
K = K{k) the complete elliptic integral of the first kind, and the modulus k depends on the 
aspect ratio of A and is here solution of K{k)/K{y/1 — k"^) — i. Using the mapping (H), the 
local rescaling factor in equation (||) is obtained, w'{z) ~ ^ [(1 — z^)(l — /c^z^)]~^/^, and inside 
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the square, keeping wi fixed, the two-point correlation function becomes (see e.g. Rcf. |2^ ) 

• 4)sq. - [^z]./^\l-z^\.\l-k^z^\) (7) 
^ ^ ^ 

k{w) 

with z{w) given by equation (^. This expression is correct up to a constant ampHtude 
determined by k{wi) which is kept fixed, but the function ip{u!) is still varying with the 
location of the second point, w. In order to cancel the role of the unknown scaling function, it 
is more convenient to work with a density profile m(w) in the presence of symmetry breaking 
surface fields hg^ on the boundary dA of the lattice A. This one-point correlator scales in the 
half-infinite geometry to(z) = {(^^ ■ ^aA(z))ioo ~ const xy^^'^ and it maps onto 

m{w) = {(fni ■ hoA{w))sq. = const x[k(w)]~'''' (8) 

where the function k{w) defined in equation (|^) again comes from the mapping. 

The application of the simple power-law in equation (|^) requires a relatively precise numer- 
ical determination of the order parameter profile of the 2d XY model confined inside a square 
with fixed boundary conditions playing the role of ordering surface fields /iaA(«,). In practice, 
no magnetic field is applied, but the symmetry is broken by keeping the boundary spins of the 
square fixed, e.g. (1)^ = (1,0), Vw G dA during the Monte Carlo simulations (in principle no 
update of these spins is performed). In the low-temperature phase, local Metropolis updates 
of single spins are known to suffer from the critical slowing down, since the system contains 
arbitrarily large clusters in which the spin orientations are strongly correlated. Statistically 
independent configurations can be obtained by local iteration rules only after a long dynamical 
evolution which needs a huge number of MC steps, sometimes beyond nowadays computers 
capability. The resort to cluster update algorithms (like Wolff algorithm) is more convenient. 
The central idea of cluster algorithms is the identification of clusters of sites using a bond 
percolation process connected to the spin configuration. The spins of the clusters are then 
independently updated. A cluster algorithm is particularly efficient if the percolation threshold 
coincides with the transition point of the spin model, which guarantees that clusters of all sizes 
will be updated in a single MC sweep. The percolation process associated to Ising or Potts 
models is the random graph model of the Fortuin-Kasteleyn representation and its threshold 
is known to coincide with the critical point of the spin model. For the XY model, the bonds 
are introduced in the Wolff algorithm through the definition of Ising variables defined by the 
sign of the projection of the spin variables along some random direction. The percolation 
threshold for these bonds coincides with the Kosterlitz-Thouless point |23j, which guarantees 
the efficiency of the Wolff cluster updating scheme at Ikt- In the low-temperature phase, 
this algorithm could be less efficient, but nevertheless preferable to a local updating. When 
one uses particular boundary conditions, with fixed spins along some surface as required here, 
the Wolff algorithm should become less efficient, since close to criticality the unique cluster 
will often reach the boundary and no update is made in this case. To prevent this, we use 
the symmetry properties of the Hamiltonian (|l|). Even when the cluster reaches the fixed 
boundaries dA{w), it is updated, and the order parameter profile is then measured with 
respect to the common new direction of the boundary spins, m(w) — {(pw • 4'dA)sq.- The 
new configuration reached would thus correspond - after a global rotation of all the spins of 
the system to re-align the boundary spins in the (1, 0) direction - to a a new configuration of 
equal total energy and thus the same statistical weight as the one actually produced. This 
technique eventually leads to a new update at each iteration: when the cluster does not touch 
the boundaries, all the spins belonging to the cluster are updated, and in the other case the 
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situation becomes equivalent to an update of all the spins but those inside the cluster. At a 
percolation threshold, interior and exterior of a percolating cluster are both large and intricated 
and the algorithm has short time correlations. After averaging over the 'production sweeps', 
one gets a characteristic smooth profile. 
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Fig. 1. - Monte Carlo simulations of the 2d XY model inside a square 100 x 100 spins (10'' MCS/spin 
after cancellation of 10^ for thermalization, cluster update algorithm). Top: on the left, the figure 
shows the rescaled order parameter density against the rescaled distance on a log-log scale at different 
temperatures below the Kosterlitz-Thouless transition temperature and slightly above. On the right, 
we show a plot of the per d.o.f. of the power-law fits as a function of the temperature for L = 48 and 
100. The sharp change of behaviour is a signature of the location of the Kosterlitz-Thouless transition, 
above which the functional expression (g) of the density profile is no longer correct. Bottom: plot of the 
density profile m{'w) with symmetry breaking surface field and of the rescaled quantity m{w) x 
vs on a log-log scale (exactly at the Kosterlitz-Thouless transition temperature). The rescaled 

profile remains constant with a very good accuracy, indicating the absence of logarithmic correction 
for this quantity. 



A log-log plot of m(w) with respect to the reduced variable k{w) is shown in figure ^ at 
two different temperatures below Tkt — 0.893(1) (the value of Tkt is taken according to 
reference |14), roughly at Tkt, and one temperature slightly above. The first message of 
the plot is the confirmation of the functional form of equation (^. One indeed observes a 
very good data collapse of the points onto a single power-law master curve (a straight line 
on this scale). The next information is the rough confirmation of the value of the critical 
temperature, since above Tkt, the master curve is no longer a straight line, indicating that 
the corresponding decay in the half-infinite geometry differs from a power-law as it should in 
the high-temperature phase. One should nevertheless mention that this is not a technique 
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Fig. 2. - Correlation decay exponent of the 2d XY model as a function of the temperature. The insert 
shows a comparison with other results (in the vicinity of the KT point) found in the literature. 



adapted to a precise determination of Tkt, since at T = 0.9 for example, the master curve is 
hardly distinguishable from a straight hne. One can for example compute the P^r d.o.f. 
as a function of temperature. It has a very small value, indicating the high quality of the fit, 
in the low-temperature phase and then increases significantly above Tkt when the behaviour 
of the density profile is no longer algebraic. The change of behaviour in the curve gives an 
approximate location of the Kosterlitz-Thouless transition temperature, and the larger the 
system size, the better this estimate. This is illustrated in figure |l|. 

It is tempting to study the role of logarithmic corrections exactly at the Kosterlitz-Thouless 
point. For that purpose, we produce data at a temperature T = Tkt and, assuming that 
logarithmic corrections should affect the density profile with symmetry breaking surface field 
in the half-infinite geometry, we make a plot of m{w) x [k(w)J^/^ vs k^w). Surprisingly, 
this leads very accurately to a constant, as shown in figure El This is a clear evidence 
that the order parameter profile with fixed boundary conditions displays a pure algebraic 
decay at the Kosterlitz-Thouless transition point. The decay exponent, if not fixed, but 
determined numerically, leads to a quite good result (for example for a size L = 100, we get 
MTkt) = 0.250(28)). 

It is thus possible to obtain the scaling dimension Xa- (T) or the two-point correlation function 
exponent rjc{T) — 2xa-{T), in the whole critical region T < Tkt as shown in figure ^. The 
most interesting result obtained here is perhaps the fact that the values of the exponent j] 
are very easy to obtain, in spite of the reputed difficulty of this problem from MC studies. 
For example the results obtained with a system of size as small as 48 x 48 are already very 
satisfying, and to get 10 values of T (10^ MCS/spin for thermalization +10^ MCS/spin for 
computation) corresponds to a quite reasonable amount of 1 1 hours of CPU time on a standard 
PC with 733 MHz processor (the algorithm is more efficient at Tkt than below) . For the 
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largest size considered here, L = 200, a simulation at the Kosterlitz-Thoulcss temperature 
needs 18 hours (10^ + 10^ MCS/spin) which is also very fast compared to today's standard 
extensive simulations where the natural time unit is the CPU-year. We mention here that we 
do not reach the accuracy of the very precise numerical determination of the magnetic scaling 
dimension (using transfer matrix techniques) by Blote and Nienhuis | |2^ . 

Once the values of these scaling dimensions are known, we can deduce numerically the 
temperature-behaviour of the compactification parameter p, and thus perhaps conjecture the 
numerical values of other scaling dimensions associated to different operators. Up to the 
authors knowledge, such an approach has not yet been done and it could be promising for 
a better understanding of the critical phase of the Kosterlitz-Thouless transitions. From the 
numerical data, ?7cr(T) is easily expanded close to the KT point T — > Tj^r^, showing a leading 
square root behaviour ||], 

r]<y^\~ 0.231(rKT - Tf'^ ~ 0.050(rKT - Tf'^, (9) 

shown as a fit in figure ^ and from which it is easy to write the first terms of an expansion 
of the compactification radius, p ^ 2[1 + 0.462(rKT - Tf''^ + 0.320(rKT - T) + . . .] when 
T — !■ Tj^rp. This leads in principle to an approximate value of the other scaling dimensions 
through equation ^. 

BB would like to thank Christophe Chatelain for stimulating discussions and Dragi Karevski 
for a critical reading of the manuscript. 
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